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Abstract 

Sterile neutrinos with masses in the keV range are considered to be a viable candidate for 
warm dark matter. The rate of their production through active-sterile neutrino transitions 
peaks, however, at temperatures of the order of the QCD scale, which makes it difficult to 
estimate their relic abundance quantitatively, even if the mass of the sterile neutrino and its 
mixing angle were known. We derive here a relation, valid to all orders in the strong coupling 
constant, which expresses the production rate in terms of the spectral function associated 
with active neutrinos. The latter can in turn be expressed as a certain convolution of the 
spectral functions related to various mesonic current-current correlation functions, which 
are being actively studied in other physics contexts. In the naive weak coupling limit, the 
appropriate Boltzmann equations can be derived from our general formulae. 
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1. Introduction 

The problem of explaining the nature of Dark Matter is a central one for cosmology. The 
most popular attempt is to postulate the existence of a relatively heavy Cold Dark Matter 
(CDM) particle related, perhaps, to softly broken supersymmetry. However, other particle 
physics candidates can also be envisaged, and there has been a recent revival particularly in 
Warm Dark Matter (WDM) scenarios. 

While WDM is just an alternative to CDM from the Dark Matter point of view, the issue 
becomes quite intriguing when other physics considerations are added to the picture. In 
particular, WDM could be realised through the existence of right-handed sterile neutrinos 1 
(see also Ref. 2 ). This possibility may then lead to some astrophysical applications [3]. 
Moreover, if there are three sterile neutrinos in total, one for each known generation, then 
they can be combined to a minimal theoretical framework, dubbed the z^MSM in Ref. and 
used to explain also the known properties of neutrino oscillations jlj and baryogenesis E] ■ 
A phenomenologically successful implementation can be obtained provided all three right- 
handed neutrinos have masses smaller than the electroweak scale. The role of WDM is played 
by the lightest right-handed neutrino, whereas the two other ones should have masses in the 
range 0(1 — 20) GeV and be very degenerate to produce the observed baryon asymmetry 
In addition, an extension of the i/MSM by a real scalar field, inflaton, allows to incorporate 
inflation [7]. 

Apart from its mass, M s , the WDM neutrino is also characterised by its mixing angle with 
active neutrinos, 9. At present the strongest observational constraints on M s and 9 come from 
two sides: structure formation in the form of Lyman-a forest observations imposes a stringent 
lower limit on M s [HI El , while X-ray constraints exclude the region of having simultaneously 
a "large" M s and 9 |10|-|14|. More precisely, if the average momentum of sterile neutrinos at 
temperatures of a few MeV coincides with that of active neutrinos then, according to Ref. 
the WDM neutrino cannot be lighter than M s ~ 14 keV. For masses in this range the mixing 
angle cannot exceed 9 ~ 2.9 x 1(T 3 (1 keV/M s ) 5/2 [T3] . 

In this corner of the parameter space the interactions of sterile neutrinos with the particles 
of the Minimal Standard Model (MSM) are so weak that the sterile neutrinos cannot equi- 
librate in the Early Universe via active-sterile transitions pQ. Therefore, information about 
initial conditions does not get lost; initial conditions do in general play a role for the current 
abundance. Evidently, it would be convenient to fix the initial conditions at some tempera- 
ture low enough such that only the z/MSM degrees of freedom play a role. Apart from initial 
conditions, one also has to fix the values of all nearly conserved quantum numbers in the 
MSM, such as the lepton and baryon asymmetries. Alternatively, one can specify the physics 
beyond the i/MSM and thus determine the initial conditions dynamically; an example of such 
a computation can be found in Ref. where the main source of sterile neutrino production 
was associated with inflation. 

One can imagine circumstances, however, where the initial conditions only play a sub- 
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dominant role. In this case the production mechanism of the sterile neutrino WDM can be 
attributed to active-sterile neutrino mixing, as in the Dodelson-Widrow scenario 1 . The 
requirements for this situation can be formulated as follows: suppose that the initial condi- 
tion is that the concentration of sterile neutrinos is zero (which is possible if the interactions 
of sterile neutrinos with all particles beyond the z^MSM such as the inflaton are negligible) 
and that there are no significant lepton asymmetries (lepton asymmetries corresponding to 
a chemical potential ///T>10 -3 would play an important role in sterile neutrino produc- 
tion Suppose also that the two heavier sterile neutrinos, present in the z^MSM, are 
heavy enough so that their decays in the Early Universe do not produce any entropy |16| . 
Then the WDM abundance can be expressed as a function of the mass M s and the angle 
6 [Tl I1UI IT71 118j . Matching the observed abundance one gets a relation between M s and 9, 
which can be confronted with the observational bounds mentioned above. 

Approximate M s -8 relations derived along these lines have been presented in Refs. I10| 
IT71 IT8] . According to the most recent analysis JH] , 



where a dark matter abundance f^DM ~ 0-22 and a QCD crossover transition temperature 
^QCD ~ 170 MeV have been inserted. If true, the combination of the Lyman-a bounds [5] and 
X-ray bounds ^Sj mentioned above rules out the Dodelson-Widrow scenario. However, other 
production mechanisms such as resonant production due to large lepton asymmetries ^H] 
(see also Ref. |19p or due to inflaton interactions [7j are feasible |16j . 

Due to the importance of the problem our aim here is to reanalyse the M s -6 relation within 
the Dodelson-Widrow scenario. In fact, a computation of the sterile neutrino production 
rate represents a very non-trivial theoretical challenge. The reason is that the region of 
temperatures at which sterile neutrinos are produced most intensely is 



At higher temperatures their production is suppressed because of medium effects |2UI . The 
temperature in Eq. (|1.2|) is very close to the pseudocritical temperature of the QCD crossover. 
Therefore, neither the dilute hadronic gas approximation nor the weakly interacting quark- 
gluon plasma picture is expected to provide for an accurate description. 

The presence of strongly interacting hadrons at these temperatures leads to two sources of 
uncertainties. A well-known one is related to the hadronic equation of state, needed for the 
time-temperature relation in the expanding Universe, entering the sterile neutrino production 
equation. Unfortunately, experiments with heavy ion collisions cannot directly measure the 
equation of state of hadronic matter. In addition, present lattice simulations with light 
dynamical quarks involve uncontrolled systematic uncertainties, such as the absence of a 
continuum limit extrapolation. In other words, the equation of state of QCD is only known 






(1.2) 
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approximately in this temperature range and contains significant systematic uncertainties 
(for a recent discussion, see Ref. (U). 

At the same time, the equation of state does matter in the computation of the sterile 
neutrino abundance. As has been mentioned already in Ref. :1] and elucidated further in 
Refs. |17j . even the purely leptonic contribution to the abundance depends significantly on 
the effective number of massless degrees of freedom, g*, which in turn changes dramatically, 
from about 60 down to about 20, when the quark-gluon plasma cools to a hadronic gas. 
However, neither the uncertainties of the equation of state, nor the subsequent uncertainties 
in the M s -9 relation, have been exhaustively investigated in these works. 

There is also a second type of a hadronic uncertainty, which has a dynamical character. 
Sterile neutrinos can be produced in reactions of two types, the first containing leptons only 
and the second having hadrons in the initial state. The hadronic reactions were omitted in 
Refs. Processes with quarks were mentioned in Ref. JJ] (not in Ref. [IHj), but without 

an explanation of how they were treated in the QCD crossover region. 

It is this second uncertainty that is the focus of the present paper. Our goal is to set up a 
general formalism for attacking it. In a later work, a numerical analysis of the M s -9 relation 
and its uncertainties will be presented. More concretely, the current number density of the 
WDM neutrinos in the Dodelson-Widrow scenario is given by 9 2 F + O(0 A ). In this paper we 
derive an expression which allows to relate the coefficient F to a certain equilibrium Green's 
function, the so-called active neutrino spectral function, defined within the MSM, while all 
dependence on the parameters of the z^MSM appears as the prefactor 6 2 . 

It is appropriate to stress that analogous relations exist in other contexts as well. For 
instance, it can be shown that the photon spectral function computed within the MSM de- 
termines the production rate of on-shell photons and dilepton pairs from strongly interacting 
systems such as colliding heavy nuclei [22], and the production rate of active neutrino pairs 
from hot / dense astrophysical environments such as the cores of neutron stars • Neverthe- 
less, we are not aware of the existence of such relations in the present context, so we want to 
discuss their derivation in a hopefully pedagogic manner. 

Given the Green's function, it should still be evaluated. As we have already mentioned, 
this turns out to be a very difficult task, since the sterile neutrino production rate peaks at 
temperatures of the order of the QCD scale. In this temperature range strong interactions 
play a dominant role, and perturbative methods fail. In the second part of our paper, we 
thus show how the active neutrino spectral function can be related to various vector and 
axial- vector current-current correlation functions defined within high temperature QCD. Such 
objects have previously been studied with a variety of methods, such as chiral perturbation 
theory, QCD sum rules, lattice QCD, and resummed weak-coupling perturbation theory, and 
also possess independent physics applications, particularly in connection with the photon and 
dilepton pair production rate computations mentioned above. 

The plan of this paper is the following. In Sec.|2]we derive the expression alluded to above, 
expressing the sterile neutrino production rate in terms of the spectral function of active 
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neutrinos, computable within the MSM. In Sec.|3]we relate the hadronic contribution to the 
active neutrino spectral function to mesonic current-current correlation functions, which can 
be defined within QCD. We also show how the result reduces to certain Boltzmann equations 
in the naive (unresummed) weak-coupling limit. We conclude and outline future prospects 
in Sec. EJ The three appendices contain certain basic definitions for the various bosonic and 
fermionic Green's functions that appear in our study, and an alternative derivation for Sec. [21 

2. General formula for the sterile neutrino production rate 
2.1. Notation 

It is a matter of convention whether the right-handed neutrinos are represented as Weyl, 
right-handed Dirac, or Majorana fermions. Choosing here the last option, the Minkowskian 
Lagrangian of i/MSM can be written as 

C = - Nnfi Nj - - MjiV/TV/ - F aI L a 4> ajiNi - F*jNi^a L L a + £ MS M , (2.1) 

where iVj are Majorana spinors, repeated indices are summed over, Mi are Majorana masses 
that we have chosen to be real in this basis, L a are the weak interaction eigenstates of the 
active lepton doublets, F a j are elements of a 3x3 complex Yukawa matrix, <ft = iT2<p* is the 
conjugate Higgs doublet, and = (1 — 7s)/2, or = (1 +^)/2 are chiral projectors. 
After electroweak symmetry breaking, (<p) = (v/\^2,0) , we define the matrix 

(M D ) aI = ^ , (2.2) 

where v ~ 246 GeV. The various neutrino fields, active and sterile, then couple to each other, 
and diagonalising the mass matrix we can define mass eigenstates in the usual way. However, 
a sterile neutrino of type / couples to an active neutrino mass eigenstate of type a with a 
very small angle only, 

01a - E i - M ^U aa , (2.3) 

a=l 

where U aa is the mixing matrix between the active neutrino interaction and mass bases. As 
mentioned above, the phenomenologically relevant part of the parameter space corresponds 
to values 9 <C 1- Therefore, the sterile neutrino interaction eigenstate of type / is to a very 
good approximation also a sterile neutrino mass eigenstate, and we can ignore the distinction 
in the following. To fix the conventions, 1 = 1 corresponds to the lightest sterile neutrino, 
contributing to warm dark matter. 

The Green's functions that we will need involve a lot of sign and other conventions whose 
definitions are unfortunately not unique in the literature. We therefore explicitly state our 
conventions in Appendices A and B, for Green's functions made out of bosonic and fermionic 
operators, respectively. Our metric convention is (H ). 
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2.2. Derivation of the master equation 

According to our assumption, the concentration of sterile neutrinos was zero at very high 
temperatures, T> 1 GeV. Moreover, because of the smallness of its Yukawa coupling, the 
lightest sterile neutrino never equilibrated. In this section we show that these two facts allow 
us to express the production rate of this neutrino through a certain well-defined equilibrium 
Green's function within the MSM. The consideration below is very general and uses only 
the basic principles of thermodynamics and quantum field theory. In particular, it does not 
require any solution of kinetic equations, nor a discussion of coherence or its loss due to 
collisions, or the like. 

The general way we proceed with the derivation is equivalent to how fluctuation-dissipation 
relations, or linear response formulae, are usually derived (see, e.g., Refs. |24l I25j ) . There 
exists, however, also an alternative derivation, which makes more direct contact with particle 
states and the related transition matrix elements and which is also somewhat shorter. The 
price is that this derivation appears to be slightly less rigorous. Nevertheless, the end result 
is identical, so we present the alternative derivation in Appendix C. 

We disregard first the Universe expansion, which can be added later on (cf. Eq. Q2.22JI ). 
Let p be the density matrix for z/MSM, incorporating all degrees of freedom, and H the 
corresponding full Hamiltonian operator. Then the equation for the density matrix is 

^ -[*«*>]. P.4) 



We now split H in the form 



H = Hmsm + H s + Hmt j (2.5) 



where JETmSM is the complete Hamiltonian of the MSM, Hg is the free Hamiltonian of sterile 
neutrinos, and flint, which is proportional to the sterile neutrino Yukawa couplings, contains 
the interactions between sterile neutrinos and the particles of the MSM. To find the concentra- 
tion of sterile neutrinos, one has to solve Eq. Q2.4|) with some initial condition. Following pQ, 
we will assume that the initial concentration of sterile neutrinos is zero, that is 

p(0)=PMSM®|0)(0| , (2.6) 

where /5msm = ^msm ex P( — /^Amsm), = 1/T, is the equilibrium MSM density matrix at a 
temperature T, and |0) is the vacuum state for sterile neutrinos. The physical meaning of 
Eq. (|2.6|) is clear: it describes a system with no sterile neutrinos, while all MSM particles are 
in thermal equilibrium. 

Considering now Hq = Amsm + Hs as a "free" Hamiltonian, and flint as an interaction 
term, one can derive an equation for the density matrix in the interaction picture, p\ = 
ex.p(iHot) p exp(-iHot) , in the standard way: 

i^ = [H 1 (t),p l (t)]. (2.7) 
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Here, as usual, H\ = exp(iHot)H- m t exp(-iHot) is the interaction Hamiltonian in the interac- 
tion picture. Now, perturbation theory with respect to Hi can be used to compute the time 
evolution of pi; the first two terms read 

pi(t) =po-i tdt' [Hi(t'),p ] + (-if tdt' f dt" [Hi(t% [Hi(t"),p ]} + ... , (2.8) 
Jo Jo Jo 

where po = p(0) = pi(0). Note that perturbation theory with Hi breaks down at a certain time 
t ~ t cq due to so-called secular terms. After t eq sterile neutrinos enter thermal equilibrium and 
their concentration needs to be computed by other means. For us t <C i e q and perturbation 
theory works well. 

We are interested in the distribution function of the sterile neutrinos. It is associated with 
the operator 

dNj _ 1 f „ 

T/ a /;q,s a /;q,s ' 



d 3 xd 3 q V 



(2.9) 



s=±l 



where a\ s is the creation operator of a sterile neutrino of type /, momentum q, and spin 
state s, normalised as 

{a /;PjS ,4 ;qt } = ^ 3 )(p - qjSuSst , (2.10) 

and V is the volume of the system. Then the distribution function dA / 7-/d 3 xd 3 q (number of 
sterile neutrinos of type / per d 3 xd 3 q) is given by 



diVj(a;,q) 
d 3 x d 3 q 



TV 



dN, 



d 3 xd 3 q 



(2.11) 



One can easily see that the first term in Eq. (|2.8f) does not contribute in Eq. ()2.11j) since Hi is 
linear in a\ s and a 7 . q s . Thus, we get that to O(0 2 ) the rate of sterile neutrino production 
reads 

= E 4 ;q ,A; q , s //t , [Hl(t) ! [ J ffl(t , ),Po]]} ■ (2.12) 

For small temperatures T <C Myy, the Higgs field in Hi can safely be replaced through its 
vacuum expectation value, so that Eqs. (|2.1|) . 1)2. 2j) imply 



Hi 



d 3 x 



{M D ) aI § a a R N I + {M* D ) aI N iaL u 



(2.13) 



where now Nj is a Majorana spinor field operator. The Nj can be treated as free on-shell 
field operators and can hence be written as 



/ = = E [a I . p , s u(I ]P ,s)e- lP - x + a\ v(I-p,s)e lP - x 
■ /(27r) 3 2p° ,=±i 

d 3 i 



/ / * n E [4 ;P ,.«( J ;p» s ) e<i '"* + a /;p,.«( I ;p. 

J ■ / (27r) 3 2p u s =±i 



s e 



(2.14) 
(2.15) 
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where we assumed the normalization in Eq. (|2.1U|) . and p° = Ep^ = p 2 + Mf , P = 
(p°, p). The spinors u, v satisfy the completeness relations J2s U (.I'> P> s )^(^j P> s ) = i> + Mi, 
Y^, s v(I; p, s)v(I; p, s) = j> — Mj, and their Majorana character requires that u = Cv T , 
v = Cu T , where C is the charge conjugation matrix. Inserting the free field operators into 
Eq. (|2.13|) . we can rewrite it as 



Hi 



d J x 



d 3 



P 



(2vr) 3 2p s =±i 



I;p,s Jl;pA X ) eiP ' X + J I;p,s( X ^ ^P' s e 



ft 



-iP-x 



where 



Ji;p,s(x) = -(M D ) aI v a (x)a R v(I\Y>, s) + (M* D ) aI u(I;p, s)a L i> a (x) . 
It remains to take the following steps: 



(2.16) 



(2.17) 



(i) We insert Eq. H2.16j) into Eq. (J2.12|) and remove the sterile neutrino creation and anni- 
hilation operators, by making use of Eq. (|2.10l) . 

(ii) This leaves us with various types of two-point correlators of the active neutrino field 
operators. We now note that correlators of the type (vp{x')v a {x)) and (vp(x')C> a (x)), 
where (...) = Tr [/3msm(—)1 an( ^ we have generalised the notation so that a, /3 incorporate 
also the Dirac indices, vanish, since lepton numbers are conserved within the MSM. 

(iii) The non-vanishing two-point functions contain the spinors u, v in a form where the 
standard completeness relations mentioned above can be used. The mass terms Mj 
that are induced this way get projected out by aL,a^. 

(iv) Introducing the notation in Eqs. (|B.1|) . (|B.2|) . the remaining two-point correlator can 
be written as 



{i> a (x')v p(x) + vp(x')v a (x)) 



d 4 P 
(2vr) 4 



-iP-(x-x') 



n>«(-P)-n< (P) 



(2.18) 



There is another term with the same structure but with x <->• x'. 



(v) It remains to carry out the integrals over the space and time coordinates. Taking the 
limit t — * oo, they yield 



lim I d 6 x /d 3 x 

t— >oo 



dt' 



e i(Q-P)(x-x>) + e i(P-Q)-(x-x') 



y(27r)V 4) (P-Q) , (2.19) 



which allows to cancel 1/V from Eq. ()2.12|) and remove P-integration from Eq. (|2.18j) . 
As a result of all these steps, we obtain 

^dV = T^o^W^WTrf^a.^^-Q) -Il< p (Q)]a R } , (2.20) 
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where we have returned to the convention that a, (3 label generations, and have expressed the 
Dirac part through a trace. Inserting Eq. ()B.8|) ; making use of the fact that 1 — np{— q°) = 
np(q°), where np(x) = l/[exp(/?x) + 1]; and observing that lepton generation conservation 
within the MSM restricts the indices a, (3 to be equal, we finally obtain the master relation 



diVj(x,q) 
d 4 xd 3 q 



i?(T,q) 



2n F (q°) 



(2^) 3 2g° 



} QY.\ M D\li^{Qa L [p aa {-Q)+ Paa {Q)]a R } , (2.21 



a=l 



where p is called the spectral function (Eq. (|B.3|0 . We stress again that this relation is 
valid only provided that the number density of sterile neutrinos created is much smaller than 
their equilibrium concentration, which however is always the case in the phenomenologically 
interesting part of the parameter space, at least for 1 = 1. 

In an expanding Universe, with a Hubble rate H, the physical momenta redshift as q(i) = 
q(to) a(to) / a(t) , where a(t) is the scale factor. This implies that the time derivative gets 
replaced with d/dt = d/dt — Hqid/dqi [22], and Eq. (|2.21j) becomes 

diV/(x,q) 



d__ H d 

dt 1 dqi 



d 3 xd 3 q 



R(T,q) 



2.22 



3. Hadronic contribution to the active neutrino spectral function 
3.1. Notation 



As stated by Eq. (|2.21|) . we need to estimate the active neutrino spectral function within the 
MSM. Given that higher-order corrections can be important, this task has to be consistently 
formulated within finite-temperature field theory. There are in principle two ways to go 
forward, the real-time and the imaginary-time formalisms. We follow here the latter since it 
can be set up also beyond perturbation theory. 

Within the imaginary-time formalism, the spectral function can be obtained through a cer- 
tain analytic continuation of the Euclidean active neutrino propagator. With the conventions 
specified in Appendix B, we denote the Euclidean propagator by (qo,q) (cf. Eq. (|B.7jO , 
Carrying out an analytic continuation, we define 



Il%J-i[q ±iO+],q) = Relr* (g°,q) ± ilmll^cAq) 



(3.1) 



where 11^ is called the retarded Green's function (cf. Eq. (|B.4|) ). The relation shown in 
Eq. (|3.1|) follows from the spectral representations in Eqs. (|B.9|) . ()B.11|) . Note that the 
imaginary part in Eq. ()3.1|) is defined as the discontinuity across the real axis: 



Imn^( g °,q) 



1 

2i 
1 

2i 



-m sc n^(-iq°, q ) 



n*(-i[q° + i0+], q) - U^(~i[q° - i0+], q 



(3.2) 

(3.3) 
(3.4) 
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where the last step introduced the spectral function (Eq. ()B.3jO and made use of Eq. (|B.12|) . 
As the imaginary-time neutrino propagator is time-ordered by construction (cf. Eq. (|B.7[) ). 
we can use functional integrals for its determination, whereby operator labels can be dropped 
from the fields from now on. 

In order to compute the Euclidean propagator, from which the spectral function follows 
through Eqs. (|3.1|) . (|3.4|) . we need to define the Euclidean theory. Given that we are interested 
in low temperatures, we can work within the Fermi- model. The interactions of the active 
neutrinos with the hadronic degrees of freedom on which we concentrate in this paper, are 
then contained in the Lagrangian 

C E = 2\[2G F (uaj^aLla + l a %aL^a + ^ v a %a L u a , (3.5) 

H ™ = d'/3B% a LUf3B , = UpBlfiaLd'ps , (3.6) 
trZ - ~ (I 4x w 75 \ -j _ ( 1 2x w 7 5 \ 
Hp = In ( 2 3 2 J U,3B /3B 3 ~2 ) 138 ' ' ' 

where a, f3 are generation indices, B is a colour index, x w = sin 2 6 W , and the Fermi-constant 
reads Gf = /4y / 2m 2 y . All fermions are Dirac fields. The fields d'p B are related to the 
mass eigenstates dps with the usual CKM matrix. The Euclidean 7-matrices are defined by 
70 = 7°, li = -*7 l > and satisfy 7* = 7^, {7 M ,>} = 25^. We have defined 75 = 70717273 = 
Z7°7 1 7 2 7 3 . Repeated /i-indices are summed over, and that they are both down reminds us 
of the fact that we are in the Euclidean space-time. We also denote Q = q^j^ and note 
that if we carry out the Wick-rotation % — > —iq° (cf. Eq. (|3.1j0 and simultanously decide to 
express the result in terms of Minkowskian rather than Euclidean Dirac- matrices, then 

i® =<foV- ( 3 - 8 ) 



3.2. General structure of the active neutrino spectral function 

Suppose now that we compute the full Euclidean neutrino self-energy within the MSM. Since 
only left-handed neutrinos experience interactions in the MSM, we expect the corresponding 
Euclidean action to have the structure 

S E = ^ E V«{Q)a R [iQ + i% aa{Q)}a L Va(Q) , (3.9) 



Qi 



a=l 



where we have defined the Fourier transforms as v a (x) = Jgexp(iQ • x)v a (Q), v a {x) = 
Jgexp(— iQ ■ x)i? a (Q). To keep the future expressions more compact, we have also factored 
the chiral projectors outside of S, but their existence should be kept in mind in the following. 
With the conventions of Eq. (|B.7|) , this leads to the Euclidean propagator 

tie in, 1 >® + g (g) ,, , m 

n m (Q) = a L _^ +i$( _ Q) a R = aL [0 + m)]2 a R , (3.10) 
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where we have made use of the property ^5 (— Q) = — % (Q), following from hermiticity (or, 
to be more precise, the so-called 75-hermiticity that replaces hermiticity in the Euclidean 
theory: the Dirac operator D satisfies 75-^75 = D) and CP-invariance. We have also left 
out the flavour indices from E to compactify the notation somewhat. 

Defining now, in analogy with the Wick-rotation of Q, a four-vector S M = (iEo,Xi); re- 
calling that Q 2 = q^q^ = —Q 2 ; and making use of Eq. (|3.8|) . we can write a general analytic 
continuation of U.^ a (Q) as 

n£,(«°,q) = nf a H<Aq) = ^ Q2 + ^^ } % (Q) a fl . (3.11) 

Writing finally ± iO + ,q) = Re£((7 ,q) ± iImE(g°,q) in analogy with Eq. mak- 
ing use of Eq. (|3.4j) . and employing the symmetry properties ReS(— Q) = — ReS(Q), 
ImS(- Q) = ImX(Q), the master relation of Eq. (|2.21|) becomes 

diV/foq) _ 4n F (g°) * \M D \ 2 aI 

d 4 xd 3 q (27r) 3 2g ( ^ 1 {[Q + ReS] 2 -[ImS] 2 } 2 +4{[Q + ReS]-ImS} 2 l ' J 

xTr{$a L (2[Q + Re£] -ImE[$ + Re2] - {[Q + Re £] 2 - [Im £] 2 } Im X )a R } , 

where S = S aQ (Q), and Q 2 = Mj. 

We remark that the trace over Dirac matrices on the latter row of Eq. (|3.12j) could trivially 
be carried out. Given that this does not simplify the structure in an essential way, however, 
we do not write down the corresponding formula explicitly. We also note that X can contain 
two types of Lorentz structures, X (Q) = Q fi(Q 2 , Q -u) + ft f2(Q 2 , Q • u), where u = (1, 0) is 
the plasma four- velocity (2Z1- However, we do not need to make a distinction between these 
two structures here. 

Now, in the absence of (leptonic) chemical potentials |15| . it is easy to see that ReX 
gets no contributions at 0{gl J /m 2 y), corresponding to 1-loop level within the Fermi model. 
The dominant contributions are 0{g^ v /m^), and originate from 1-loop graphs within the 
electroweak theory [281 129j . The dominant contributions within the Fermi-model are of 2- 
loop order, 0(g^ a /my ir ), and thus suppressed with respect to the 1-loop effects from the 
electroweak theory. In contrast, Im S requires on-shell particles on the inner lines, and cannot 
at low energies E <C my/ get generated within 1-loop level in the electroweak theory (more 
precisely, ImS is exponentially suppressed by ~ exp(— mw/T)). The dominant contributions 
are 0(g^ J /my V ) and can be computed within the Fermi-model. It is these contributions that 
we concentrate on in the following. 

For general orientation, it is useful to note that if we assume Re X, ImS <C Q, as is 
parametrically the case at low energies, then Eq. (J3.12|) can be simplified to 

dJV f (x,q) ^ 4n F (g°) * I^d|^ r 1 , . 

"dWq"- (2l^^ 1 ^ 2 - Tr ^ aiIm ^H ' < 3 ' 13 ) 



10 



Given that the large-time decay of the retarded propagator H R is determined by the structure 
g° + iImS° (cf. Eq. flTTTR and given our conventions in Eq. (|B.4j) . we expect the behaviour 
H R (x°) ~ f a exp(— iq°x°)/(q° + ilm£°) ~ exp(— x° ImS°). Therefore ImS° has to be 
positive; in fact, we can define ImS° = r„/2, where T v is called the active neutrino damping 
rate. As Eq. (|3.13[) shows, Im£° > also leads to a positive sterile neutrino production rate. 

3.3. Relation of active neutrino and mesonic spectral functions 

With the conventions set, we need to determine i^ aa (Q). A simple computation to second 
order in the Fermi interaction (as already mentioned, the first-order contribution vanishes in 
the absence of chemical potentials) yields 

4 aa {Q) = AG 2 F £ 4 PH 2 >qg,(£), (3-14) 
H=w,z T n b (Q + R) +m lH 

where pw = 2, pz = 1/2 are the "weights" of the charged and neutral current channels; 
mi w = mi a is the mass of the charged lepton of generation a; m\ z = m Ua = is the mass 
of the MSM active neutrino; = (fo,r) denotes bosonic Matsubara four-momenta; and 
we have defined the Euclidean charged and neutral current correlators in accordance with 
Eq. (|A~8l) . viz. 



C™{R) = f e^-y\H^(S:)HT(y)) , C%{R) = f e^M/flJ (x)H?(y)) , (3.15) 

•J X J X 

where x^ = (x°,x l ), f M = (fo,fi) = (fo, — r l ), x-R = x^f^ = x°fo—x l r l , and /- = J^dx Jd 3 x. 
The tildes in C's and -fTs remind us of the fact that we are for the moment using Euclidean 
Dirac-matrices in the hadronic currents. We also point out that the Dirac algebra remaining in 
Eq. (|3.14|) cannot in general be greatly simplified, since the functions C^ U {R) can in principle 
contain both symmetric (e.g. R^Ru) and antisymmetric (eapfiv'UaRp) structures in fi <-» v. 

It is important to stress now that in Eq. Q3.14|) we have assumed a free form for the 
lepton/neutrino propagators inside the loop. Naturally, one could also allow for a general 
structure, such as the one in Eq. (|3.10|) . to appear here. To keep the discussion as simple as 
possible, however, we treat the inner lines to zeroth order in Gf in this paper. 

To see how Eq. (|3.14[l can be analysed, let us simplify the notation somewhat. The essential 
issue is what happens with the Matsubara frequencies, and we hence rewrite the structure as 

\ _ v/ Uu(i[qo + r \) *h r~ \ Q /"a^ /qi«n 

l% (90 ' q) = % few ^ (ro ' r) ' V J ' ( } 

R b Q+R 

where 



E 1 = y /( q + rf + mf H . (3.17) 

The drawing in Eq. (|3.16j) illustrates the momentum flow, with the thick line indicating the 
composite mesonic propagator. The challenge is simply to rewrite Eq. (|3.16j) in a form where 
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the analytic continuation needed for Re £ and Im £ can be carried out in a controlled way, 
without generating any non-converging sums or integrals. 
In order to proceed, we first rewrite Eq. (|3.16[) as 



-q -r _ 2 p 2V" ' ^ 
S + £j x 



(3.18) 



where f T = Jd 3 r/(27r) 3 , and f b,sof denote bosonic and fermionic Matsubara frequencies, 
respectively. Furthermore, we note that the Kronecker <5-function can be expressed as 



5- - ~ — T A T J T ( s o-qo-r ) 



Thereby the correlation function becomes 



i% (?o>q) 



dre"^ 



r JO 



- s of 



(3.19) 



(3.20) 



The sum inside the first square brackets can be performed according to Eq. (jB.14|) . In order to 
handle the second square brackets, we express Cfj^(ro, r) through the spectral representation 
in Eq. ()A.14j) . Inserting this into Eq. (|3.2U|) and changing orders of integration, we obtain 



it (90, q) 



2E 1 



7T ^ 



dr e 



-trq 



' ob 



uj — iro 



(3.21) 



As the next step, the sum in Eq. (|3.21j) can be performed. In fact, the result can immediately 
be obtained by choosing suitable constants c,d in Eq. (|A.17|) : for < r < (3, 



-irro 



'ob 



uj — iro 



(3.22) 



where ne(x) = l/[exp(/3x) — 1]. The integral over r can then be carried out according to 
Eq. (|B.15|) . leading finally to 



(?o,q) 



n F (^i) 



oo vr 
e P(u+Eh.) _|_ 1 

ig + w + Ei 



/»2/(w,r)ra B M x 



e f3w _|_ e /3Ei 
i 

ig + w - -Ei 



(3.23) 



It remains to carry out the analytic continuation go — * — ±i0 + ] and to take the real and 
imaginary parts. In particular, taking the imaginary part as defined by Eq. (|3.3|) goes with 
Eq. (|A.12|) . given that % only appears in a simple way in Eq. (|3.23j) . Carrying then out the 
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integration over uj to remove the <5-functions, and returning simultaneously to Minkowskian 
Dirac-matrices (tuTu = 7 M 7/x)) whereby the tildes can be removed, we arrive at 



x v*, n s f cosh(/?g°/2 



r 4E 1 cosh(/3^i/2) 
f^(- El )p^(- q ° - Ek, r) f^jEQp^-q + E u r) 
sinh[/3(g° + E x )/2] sinh[/?(g° - E{)/2] 



(3.24) 



Defining 



A(p°,p,m) = f +m , (3.25) 



and reintroducing masses in the numerators to remind us of the fact that the A-functions 
are to be evaluated at the on-shell points (the masses are in any case deleted by the chiral 
projectors), we can return to the complete notation: 



T v* / o n a„2 f d3r cosh(/3<?°/2) 



H=W,Z 



(2vr) 3 4Ei cosh(/3Ei/2) 

(3.26) 



S inh[/3(,o + jBl )/ 2 ] ^ ( -« 

Here E\ is from Eq. (|3.17|) . Let us note that p^ u vanishes at zero frequency, so that the poles 
originating from the sinh-functions in Eq. ()3.26jl are harmless. 

To conclude, ImS can be expressed as a three-dimensional spatial integral, with certain 
hyperbolic weights, over the mesonic spectral functions related to the charged and neutral 
currents. On the other hand, Re E requires a four-dimensional principal value integration over 
the same spectral functions, as dictated by Eqs. ()3.23|) . (|A.12|) . 

3.4. Reduction of mesonic spectral functions 

If we make certain assumptions about the quark mass matrix, the mesonic correlators in 
Eq. (|3.15|) . as well as the corresponding spectral functions p^ v that appear in Eq. (|3.26|) . can 
be reduced to a small set of quantities that have been widely addressed in the literature. 

To start with, it is probably a good approximation at temperatures 100 MeV < T < 400 
MeV to treat the three lightest quarks as degenerate, with a certain mass m q , while the three 
heaviest quarks can be assumed infinitely heavy, and ignored. In this limit the theory has 
an exact SUy(3) flavour symmetry, which guarantees that we can split the correlators into 
flavour singlets and non-singlets. Defining T a to be traceless and Hermitean, and T° to be 
the 3x3 unit matrix, we can then construct the flavour non-singlet and singlet vector and 
axial currents, 

V? = $%T a ^, K° = ^TV, A a =^ 5 T a il>, A°= i;% l5 T i; . (3.27) 
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In general, the flavour symmetry allows for six correlation functions: 

Tr [T a T b ]C^ u (R) = J_ e^* (v« (x)V b (0)) , (3.28) 

Tr [T a T b ]C'^ u (R) = J_ e^(i» (5)i*(0)) , (3.29) 

Tr [T a T b ]C™(R) = jf e iM {v^x)A b M + A«(x)V b (0)) , (3.30) 

C™(R) = J/ A %V°(x)V u °(0)) , (3.31) 

C$(R) = j£e^(i°(£)i°(0)), (3.32) 

<5™°(£) = |_e^(y M (x)i°(0) + i°(x)K (0)). (3.33) 

However, only four among these are non-trivial in the mass-degenerate limit that we are 
considering: the vector and axial currents have opposite transformation properties in charge 
conjugation C, which implies that Cj^}(R), Cj^°(R) vanish in this case. 

Now, the correlators that appear in Eq. (|3.15|) can be expressed in terms of the ones just 
defined. We obtain 

CZ{R) = ' Kd|2 ''V VUS? 1 \CUR) + C*(R)] , (3.34) 



[IV 1 



3 



\2/W / 6\ , r<A 



CU R ) = o (l-2x w yC^(R)+C^(R) +- C™(R) + C™(R) , (3.35) 



1 



36 



where V$ are elements of the CKM matrix. Identical relations hold for the spectral functions. 

Under further assumptions, the set of independent correlators can still be reduced. In 
particular, taking the chiral limit m q — > 0, the Ward identity following from applying an 
infinitesimal non-singlet axial transformation on the correlator (V^(x)A b (0)} states that 
CY,(R) = Cp V (R). For the singlets this is not true in general, in spite of the fact that 
anomalous processes are in some sense suppressed at high temperatures (see, e.g., Ref. |50*|). 
So in principle there remain three independent functions to determine, CY,, CY®, Cm?. 

Now, in the three-flavour theory, the non-singlet vector correlator C^ v determines directly 
the photon spectral function [22], and is thus relevant for computing the photon and dilepton 
production rates, which are among the prime observables for heavy ion collision experiments. 
Therefore it has been addressed with a variety of methods in the literature, starting with 
various perturbative treatments [22] as well as with so-called thermal sum rules [HJ- The 
(resummed jH2]) perturbative treatments have reached a great degree of sophistication by 
now [HS], with different strategies applicable in different parts of the phase space. On the 
other hand, it has also been realised that an analytic continuation of imaginary time cor- 
relators defined on the finite r-interval can be carried out in principle |34j . which opens up 
the possibility of lattice QCD determinations. There have indeed been attempts at practical 
implementations of a certain analytic continuation from numerical data [33], though they 
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are not without problems for the moment |36| . Finally, for T <C 150 MeV, chiral pertur- 
bation theory can be systematically applied jSZ] or at least used as a solid baseline for the 
computation of C^ u (and thus of p^ u ) [38] . 

The vector singlet C*^ , on the other hand, can be associated with the baryon num- 
ber current, whose susceptibility \ (which is an integral over the spectral function, x = 
J^^dio Pqq (uj , 0) / ttlu) is argued to be relevant for so-called event-to-event fluctutations in 
heavy ion collision experiments [30] • In resummed perturbative treatments 0Hj the differ- 
ence between the vector singlet and non-singlet susceptibilites (this difference is often called 
the "off-diagonal quark number susceptibility") is however very small, being suppressed by 
af ln(l/a s ) [3Tj, so that at high enough temperatures C^i? can well be approximated (up to 
an overall factor) by C^. At lower temperatures close to T ~ 150 MeV, on the other hand, 
a lattice determination would again be needed; unfortunately, the difference between C^® 
and is given by a disconnected quark-line contraction which is technically rather difficult 
to measure accurately |15) . We should of course stress that the susceptibility alone contains 
much less information than the complete function CY£, or its analytic continuation pff. The 
general high-temperature structure of the latter has been analysed in Ref. [33] • Finally, at 
very low temperatures, chiral perturbation theory predicts that is strongly suppressed 
with respect to C^ u . 

Much the same comments can be made for the axial singlet current, Cff. At high enough 
temperatures, where resummed QCD perturbation theory is applicable, it agrees up to a 
certain order in the resummed perturbative expansion with CY^ . At lower temperatures, 
the difference between Cffi and becomes significant. This difference is of course quite 
interesting in its own right, being related to the r/-meson and the chiral anomaly. However, 
analytic and numerical treatments are demanding in this range. At very low temperatures, 
chiral perturbation theory indicates that Off remains more significant than CY£. 

3.5. Perturbative limit 

Returning finally to the simplest possible logic, we wish to inspect the mesonic spectral 
function p^ u in naive perturbation theory. The purpose is to show that in this case, Eq. (j3.2(i|) 
leads to the familiar structure of Boltzmann equations. (The procedure we follow is analogous 
to the one first worked out for simpler cases in Ref. |44j. and also used in previous active 
neutrino damping rate computations |45j.) 

In order to reach this goal, we write generic hadronic currents (Eqs. (|3.6|) . Q3.7|) ) in the 
form Hfj,(x) = ^{x)^ q^{x) , where q2, denote quark fields and T is some Dirac matrix 
structure. Carrying then out the contractions in the correlators of Eq. (|3.15j) . we obtain 

^ H ~ AT ^ f -i(f +ff)+m 2 . -if +m 3 . ) f . 
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where N r = 3 is the number of colours, and 



E 2 = yj(t + r) 2 + m\ , E 3 = ^t 2 + m 2 . (3.37) 

The important issue is again what happens with the Matsubara frequencies. Following the 
steps in Section 13.31 and denoting the complete numerator of Eq. ()3.36|) by a function g^ u , 
we obtain 



9^{ito + if ,it ) 

••nf 



UO—tQ—fo r~2 



*0f'"0f 

/9 



[ul + El][t> + El 



[u 2 + £ 2 2 ^ 2 + £ ; 



"(if 



of 



n F (E 2 )n F (E 3 
AE 2 E 3 



dr e 



" irfo +^(-^2,+^ 3 )e (/3 - r)(£;2+£;3) - 

-~g^(+E2,+E 3 )e^ E * +TE > + 
+ g^(+E 2 ,-E 3 )e^+ E ^] , 



(3.38) 



where we made use of Eq. (jB.14|) . The integral remaining can be carried out by using 
Eq. (|A.18|) (recalling that r"o is bosonic), and the spectral function follows then by picking up 
the discontinuity across the real axis, according to Eq. (|A.12|) : 



-7T 



t AE 2 E 3 



+5{r° + E 2 + E 3 )g^{-E 2 , +E 3 )(1 - n F2 - n F3 ) + 

+ S(r° + E 2 - E 3 )g^(-E 2 , -E 3 )(n F2 - n F3 ) + 
+ 5(r° -E 2 + E 3 )g^(+E 2 , +E 3 )(n F3 - n F2 ) + 

+ 5(r° -E 2 - E 3 )g^(+E 2 , -E 3 ){n F2 + n F3 - 1) 



(3.39) 



where we have denoted nFi = n F (Ei). 

The expression in Eq. (|3.39|) can now be inserted into Eq. ()3.26j) . There are eight different 
terms in total. The hyberbolic functions in Eq. (|3.26f) can be rewritten in terms of n-p's and 
ne's, and making use of relations such as 

S(-q° +E X +E 2 + E 3 )nv{q® - E x )(l - n F2 - n F3 ) = 5(-q° + E X +E 2 + E 3 )n F2 n F3 , (3.40) 

they can be reorganised in a simpler form. In view of Eq. H2.21j) . we also choose to factor 
out rip 1 (q ) from all the terms. After a number of trivial if tedious manipulations, we finally 
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arrive at 



d 3 i 



d 3 



P.3 



H=W,Z 



(2vr)V 4 

+ (2vr) 4 5^ 4 

+ (2^)V 4 

+ (2vr) 4 5^ 4 

+ (2^)V 4 

+ (2^)V 4 

+ (2^)V 4 

+ (2^)V 4 



Pi+ P2+P3- Q) n F1 n F2 n F3 A(-mi H ,m 2 , -m 3 ) + 2 -A-*. . 

3' 

-P2 + -P3 - -Pi - Q) n F2 n F %(l -n F i)i(m lff ,m 2 ,-m 3 ) + 
P1 + P3-P2- Q)n F in F3 (l - n F2 )^(-m; H , -m 2 , -m 3 ) + \£ 
Pi + P 2 - P3 - Q)ra F in F2 (l - n F3 )^(-m^,m2,m 3 ) + ^X* 
Pi - P2 - P3 - Q)n F i(l - n F2 )(l - n F3 )^l(-m iH , -m 2 ,m 3 ) + 1 



P 2 - Pi -P3 - Q)n F2 (l - n F i)(l - n F3 )„4(m; H ,m 2 ,m 3 ) + 
P3 - Pl -P2 - Q)n F3 (l - n F i)(l - n F2 )A(mi H ,-m 2 , -m 3 ) + 3 
-Pi - P2 - P3 - Q) (1 - ifi)(1 ~ n F2)(l - n F3 )_4(m; H , -m 2 ,m 3 ) 



2 

- Q 
3 



where 



A{m lH , m 2 , m 3 ) = 7 M 0f 1 + "^h* Tr (f 2 + m 2 )j tl T(f 3 + m 3 ) 7l T 



(3.41) 



(3.42) 



Here Pj = (Pj,pj) are on-shell four- momenta, and the energies have been transformed from 
Eqs. (|3~T7j> . (|3~37|l into 



Pl = VPl + ™/ H , p2 = VP2 + m 2> ^3-V P 3 +m 3- 



(3.43) 



The graphs in Eq. (|3.41() illustrate the various processes, with time and momenta assumed 
to run from left to right, and arrows indicating particles / antiparticles in the usual way. 
The general structure is clearly what we expect from Boltzmann equations, however we have 
arrived at it without any model assumptions, evaluations of scattering matrix elements, or 
spin averages. It is easy to check (by making use of the properties of the n F 's as well as 
the fact that the first argument of A(mi,m 2 ,m 3 ) can be dropped as it will in any case be 
projected out by ol, or, and that -4(0, — m 2 , — m 3 ) = .4(0, m 2 ,m 3 ) for T = a + 675) that the 
expression in Eq. (|3.41|) is symmetric in Q — > —Q, as it should be. 

For illustration it is useful to take one more step, and insert this result into the simplified 
form in Eq. (|3.13j) . Let us consider the contribution from the third term in Eq. ([3.41)1 . for 
instance, and choose the charged current, i.e. H = W. Then T = cll, and we obtain 



dJVj(x,q) = 16N C G 2 F ' \M D \lj f d 3 Pl 
d 4 xd 3 q (2vr) 3 2g° 



, MI 



d 3 P2 



(2vr) 3 2Pi J (2vr) 3 2P 2 J (2vr) 3 2P 3 



d 3 P3 
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(2^) 4 ^ 4 )(P! + P3-P2- Q)n F (u ■ Pi)n F (n • P 3 )[l - n F (u ■ P 2 )] x 



x Tr 



Tr 



*2 - "^h^Cf 3 ~ m 3 )j v a L 



(3.44) 



where = has only been kept to formally indicate the anti-particle direction of the line, 
and we have written everything in a Lorentz-covariant form. All explicit masses drop out 
because of the chiral projectors. The Dirac algebra is elementary, 



Tr 



Q-ffi-fan Tr f 2 l„fzlua L = 16 Q • P3P1 • P; 



(3.45) 



This expression is positive, and carrying out the phase-space integral in Eq. (|3.44[) . one 
obtains a finite function of |q|. Let us remark, though, that while the integration over the 
energy-conservation constraint is simple in the center-of-mass frame, the plasma four-velocity 
u becomes non-trivial if this frame is chosen; in general, therefore, the phase-space integrals 
that appear in Eq. I|3.44[) are technically non-trivial. 

It is appropriate to end now by pointing out that Eqs. (|3.44j) . (|3.45|) contain significant 
differences with respect to the spin-averaged matrix elements squared that appear in ac- 
tive neutrino scattering cross-sections [IS], and have sometimes also been inserted into the 
Boltzmann equations for sterile neutrino production. In particular, the left-most trace in 
Eq. Q3.45JI contains the projector an, rather than cll] this is because Eq. (|2,21l) contains the 
active neutrino propagator, rather than the self-energy. As a consequence, Eq. (|3.45|) does 
not lead to a purely s-channel momentum structure (16 QP 2 P\- P3) like the active neutrino 
scattering cross sections |45| . Whether this special example has any practical significance is 
not clear at this stage, but it illustrates the advantages of our framework where the correct 
structures are produced automatically from thermal field theory, rather than having to be 
inserted by hand. 



4. Conclusions and Outlook 

While the sterile neutrino production rate has previously been investigated in the literature 
in some detail, the hadronic contributions to it have never been analysed properly. These 
contributions involve strongly interacting dynamics at temperatures of the order of the QCD 
scale, where neither perturbation theory nor the dilute hadronic gas approximation are valid. 

To confront this situation, we have derived a general relation that expresses the sterile 
neutrino production rate in terms of the active neutrino spectral function, computable by 
using equilibrium thermal field theory within the Minimal Standard Model (Eq. (|2.21j) ). The 
active neutrino spectral function can in turn be expressed in terms of the real and imaginary 
parts of the active neutrino self-energy (Eq. I|3.12|) ). These equations show that hadronic 
contributions may play a role in three different ways: 

(1) Most importantly, the hadronic degrees of freedom contribute at leading order to the 
imaginary part of the active neutrino self-energy, ImS. Therefore, we have expressed 
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the hadronic contribution to ImE as a certain convolution of the spectral functions 
related to hadronic current-current correlation functions (Eq. 1)3.26)) ). The latter can in 
turn be expressed in terms of standard vector and axial- vector correlators (Eqs. 1)3.34)1 . 
1)3.35)1 ) that can be studied with a number of different theoretical methods, and are also 
partly related to experimental observables addressed in the heavy ion program. 

(2) The parametrically dominant contribution to the real part of the active neutrino self- 
energy, ReS, arises from 1-loop graphs within the electroweak theory, and does not 
contain hadronic effects. On the other hand, subdominant contributions, formally sup- 
pressed by a w , do contain hadronic effects. As shown by Eq. (|3.23|) . these contributions 
can be expressed in terms of a certain weighted integral over the same hadronic spectral 
functions that appear in the imaginary part. 

(3) Finally, the sterile neutrino production rate equation contains a time derivative d/dt, 
the Hubble rate H, and the temperature T (cf. Eq. ()2.22)0 . In cosmology, all of these 
are related via the Einstein equations. The relation is again sensitive to hadronic 
effects, via the equation-of-state of the primordial plasma. The current status and 
phenomenological fits for all the thermodynamic functions that appear in the time- 
temperature relation can be found in Sec. IV of Ref. |21j . 

To summarise, our formulae should allow to estimate for the first time the systematic 
hadronic uncertainties that exist in computations of sterile neutrino production through 
active-sterile transitions. A numerical evaluation of these effects is in progress. 

Apart from these non-perturbative aspects, we have also demonstrated that our general 
formulae allow to derive, without further assumptions, the appropriate Boltzmann equations 
that apply in the naive weak-coupling limit (Eq. 1)3. 41J) ) . It may in fact be useful to repeat our 
computations for the leptonic contributions as well, since a first-principles derivation frees us 
from the need to argue about spin averages or symmetry factors, and produces automatically 
the correct Dirac and chiral structures for the Boltzmann equations. 

Acknowledgements 

The work of T.A. was supported in part by the grants-in-aid from the Ministry of Education, 
Science, Sports, and Culture of Japan, Nos. 16081202 and 17340062, and that of M.S. by the 
Swiss National Science Foundation. 

Appendix A. Basic relations for bosons 

We list in this Appendix some common definitions and relations that apply to two-point 
correlation functions built out of bosonic operators; for more details see, e.g., Refs. |24~) 125]. 
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We denote Minkowskian space-time coordinates by x = (t, x l ) and momenta by Q = (q°, q % ), 
while their Euclidean counterparts are denoted by x = (r, x 1 ), Q = {qo,Qi)- Wick rotation is 
carried out by r <-> it, % <-> — Arguments of operators denote implicitely whether we are 
in Minkowskian or Euclidean space-time. In particular, Heisenberg-operators are defined as 



(A.l) 



6(t, x) = e iHt 6(0, x)e- iHt , 6(r, x) = e Hr O(0, x) e - Hr . 

The thermal ensemble is defined by the density matrix p = Z^ 1 exp(—f3H). 

We denote the operators which appear in the two-point functions by (j> a (x), (f>p(x). They 
could be elementary field operators, but they could also be composite operators consisting of 
a product of elementary field operators. 

We can now define various classes of correlation functions. The "physical" correlators are 
defined as 



h(x)4>Uo 



Kp(Q) = /did 3 xe^(^(O)0 a (x 
P*f>(Q) = |dtd 3 xe^(i[^(x),^(0) 



(A.2) 
(A.3) 
(A.4) 



where p a p is called the spectral function, while the "retarded" / "advanced" correlators can 
be defined as 



n«(Q) ee *|did 3 xe^ 

n^(Q) = ijdtd 



3 xe iQ-x 



l a (x),4>l(o)]e(t)) , 
4> a (x),4>l(o)]e(-t] 



(A.5) 
(A.6) 



On the other hand, from the computational point of view one is often faced with "time- 
ordered" correlation functions, 



n^(Q) = |dtd 3 xe^(^(x)^(O)0(t)+^(O)^(x)0(-t)), (A.7) 



which appear in time-dependent perturbation theory, or with the "Euclidean" correlator 



nf^(Q) = J q dr|d 3 xe^(4(x)^(0)), 



(A.8) 



which appears in non-perturbative formulations. Note that the Euclidean correlator is also 
time-ordered by definition, and can be computed with standard imaginary-time functional 
integrals in the Matsubara formalism. 

Now, all of the correlation functions defined can be related to each other. In particular, 
all correlators can be expressed in terms of the spectral function, which in turn can be 
determined as a certain analytic continuation of the Euclidean correlator. In order to do 
this, we may first insert sets of energy eigenstates, to obtain the Fourier-space version of the 
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so-called Kubo-Martin-Schwinger (KMS) relation: n< /3 (Q) = e'^U^iQ). Then p a p(Q) = 
IKp(Q) ~ Kp(Q)\/ 2 and > conversely, 

U> p (Q) = 2[1 + n B (q°)}p a p(Q) , Ii< p {Q) = 2n B (q°)p a( ,(Q) , (A.9) 
where n B (x) = l/[exp(/3x) — 1]. Inserting the representation 

f°° e~ iu)t 

9(t)=i — — (A.10) 

into the definitions of T1 R , IL A , we obtain 

U R (n) _ f°° duJ P^(^,q) A (n) _ f°° PapMl s 

n^^-y^V^o— o+ ' u ^ Q) -J^—u- q ° + l o+ ■ (A ' n) 

Doing the same with T1 T and making use of 

' p(1)t^(A), (A.12) 



A±i0+ VA 

produces 



n^(Q) = r^ ^^l +2p^( g °,q)n B ( g °) • (A.13) 

Finally, writing the argument inside the r-integration in Eq. (|A.8|) as a Wick rotation of the 
integrand in Eq. HA.2j) . which in turn is expressed as an inverse Fourier transform of H > (Q), 
for which Eq. (|A.9|) is inserted, and changing orders of integration, we get 



.5,(4) = \\t^ r p^u^ - r*zi*m 

'■' ' 27T P J-oo 7T lo -tqo 



oo 



This relation can formally be inverted by making use of Eq. (|A.12|) . 



1 



Pap(q ,q) = — Discn£g(9o -> -*g°,q) , (A.15) 

where the operation Disc is defined in Eq. 

We also recall that bosonic Matsubara sums can be carried out through 

nB ^) r, + d)e (^-r)E + {cE + d)e rEl j (A17) 



2£ 



where u>b = 27rTn, with n an integer, and we assumed < r < /3; and that a typical 
integration yields 

rP 1 - p-^ A 

/ dre- r ^b+A) = . ( A-lg ) 

Jo iuj h + A 
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Appendix B. Basic relations for fermions 



We list in this Appendix some common definitions and relations that apply to two-point 
correlation functions built out of fermionic operators; for more details see, e.g., Refs. [241125]. 

We denote the operators which appear in the two-point functions by v a {x), vp(x). They 
could be elementary field operators, in which case the indices a,/3 label Dirac and/or flavour 
components, but they could also be composite operators consisting of a product of elementary 
field operators. 

Like in the bosonic case, we can define various classes of correlation functions. The "phys- 
ical" correlators are now set up as 

n^(Q) = |dtd 3 xe'^(z> Q (x)^(0)), (B.l) 

Kp(Q) = fdttfxe^-tpMOM), (B.2) 

Pap{Q) = Jdtd 3 xe^ x (^{i>^x),D p (0)}), (B.3) 
where p a p is the spectral function, while retarded and advanced correlators can be defined as 

n&CQ) = iJdtd 3 xeW-*({i> a (x),0p{O)}O(t)) t (B.4) 

<j(Q) = ijdtd 3 X e iQ - x (-{u a (x),Pp(0)}e(-t)) . (B.5) 
On the other hand, the time-ordered correlation function reads 

n^(Q) = Jdtd 3 xe iQ - x (v a (x)9p(0)9(t)-9p(0)v a (x)8(-t)) , (B.6) 
while the Euclidean correlator is 

n§s(Q) = J^dr yd 3 xe^(i> a (x)^(0)) . (B.7) 

Note again that the Euclidean correlator is time-ordered by definition, and can be computed 
with standard imaginary-time functional integrals in the Matsubara formalism. 

Like in the bosonic case, all of the correlation functions defined can be expressed in terms 
of the spectral function, which in turn can be determined as a certain analytic continuation 
of the Euclidean correlator. First, inserting sets of energy eigenstates, we obtain the KMS- 
relation in Fourier-space, II^Q) = - e -^ n > /3 (Q). Then p af3 {Q) = Pjg(O) - n<,(Q)]/2 
and, conversely, 

nJs(Q) = 2[1 - n F (q°)]p a p(Q) , U<p(Q) = -2n F (q°)p a p(Q) , (B.8) 

where uf(x) = l/[exp(/?x) + 1]. Inserting the representation of Eq. (jA.lOj) into the definitions 
of Tl R , H A produces 

1 J-oo vr lo - q u - ^0 + H J-oo n u - q u + t0 + 
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Proceeding similarly with I1 T and making use of Eq. (|A.12|) , we obtain 



nL(Q) 



duj ip a p(uj,ci) 
> 7T q° — uj + i0+ 



2 Pa/3 (q°,q i )n F (q ) . 



(B.10) 



Finally, writing the argument inside the r-integration in Eq. (|B.7|) as a Wick rotation of the 
inverse Fourier transform of the left-hand side of Eq. (jB.l|) . inserting Eq. (|B,8|) . and changing 
orders of integration, we get 



^' " /i,r / n> ( W ,q 



2tt 



30 du p af3 (uj,q) 
00 vr lo — i% 



Like in the bosonic case, this relation can be inverted by making use of Eq. (|A.12[) . 



p(q°,q) = i Disc IF 5 (go 



We also recall that fermionic Matsubara sums can be carried out through 



f 2 +£ 2 



tip(S) 



W f 2 + i? 2 



2£ 



(f3-T)E 



(cE + d)e T 



(B.ll) 

(B.12) 

(B.13) 
(B.14) 



where u>f = 2i\T{n + 5 )> with n an integer, and we assumed < r < /?; and that a typical 
integration yields 

' dTe -^,+A) = l + e-^ 



/ 



zcjf + A 



(B.15) 



Appendix C. An alternative derivation of Eq. (12.211) 

We present in this Appendix an alternative derivation (following Ref. for example) of 
Eq. 1)2.21(1 . which is technically somewhat simpler than the one in the main text, but with 
the price of containing a few heuristic steps. The end result is nevertheless identical. 

The starting point is the interaction Hamiltonian in the phase with broken electroweak 
symmetry, Eq. (|2.13|) . Consider now an initial state \X) = |i) ® |0) and a final state \J-) = 
|f) (g) |I;q, s), where the right subspace contains the sterile neutrinos, and 

|I; q, s) = a}. q JO> , (/; q, s\ = <0|a J;q , s . (C.l) 

The transition matrix element can immediately be written down, 

Trx = (F\ [dtH!(t)\l) = Ut^x- ^* (f|J f;q , a (g)|i) , (C.2) 
J J V /(27r) 3 2g° 
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where we inserted Eq. (|2.16|) . and J/ ; q, s is from Eq. (|2.17[) . The production rate can then be 
obtained by summing over all initial states, with their proper Boltzmann weights, and over 
all allowed final states: 

^ s=±l f,i 



(2vr) 3 2g° 



s=±l 



where is the volume, At is the time interval, Z the partition function, we defined a thermal 
average by (...) = Z _1 Tr [exp(— (3 H^su) (■■■)): an d made use of translational invariance. 

It remains to repeat steps (ii), (iii) in the paragraph following Eq. Q2.17[l . We thus arrive 
directly at Eq. lf!OU)l . 
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